home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 22 / Cream of the Crop 22.iso / math / ast53src.zip / MATRIX.C < prev    next >
C/C++ Source or Header  |  1996-09-29  |  22KB  |  730 lines
  1. /*
  2. ** Astrolog (Version 5.30) File: matrix.c
  3. **
  4. ** IMPORTANT NOTICE: The graphics database and chart display routines
  5. ** used in this program are Copyright (C) 1991-1996 by Walter D. Pullen
  6. ** (Astara@msn.com, http://www.magitech.com/~cruiser1/astrolog.htm).
  7. ** Permission is granted to freely use and distribute these routines
  8. ** provided one doesn't sell, restrict, or profit from them in any way.
  9. ** Modification is allowed provided these notices remain with any
  10. ** altered or edited versions of the program.
  11. **
  12. ** The main planetary calculation routines used in this program have
  13. ** been Copyrighted and the core of this program is basically a
  14. ** conversion to C of the routines created by James Neely as listed in
  15. ** Michael Erlewine's 'Manual of Computer Programming for Astrologers',
  16. ** available from Matrix Software. The copyright gives us permission to
  17. ** use the routines for personal use but not to sell them or profit from
  18. ** them in any way.
  19. **
  20. ** The PostScript code within the core graphics routines are programmed
  21. ** and Copyright (C) 1992-1993 by Brian D. Willoughby
  22. ** (brianw@sounds.wa.com). Conditions are identical to those above.
  23. **
  24. ** The extended accurate ephemeris databases and formulas are from the
  25. ** calculation routines in the program "Placalc" and are programmed and
  26. ** Copyright (C) 1989,1991,1993 by Astrodienst AG and Alois Treindl
  27. ** (alois@azur.ch). The use of that source code is subject to
  28. ** regulations made by Astrodienst Zurich, and the code is not in the
  29. ** public domain. This copyright notice must not be changed or removed
  30. ** by any user of this program.
  31. **
  32. ** Initial programming 8/28,30, 9/10,13,16,20,23, 10/3,6,7, 11/7,10,21/1991.
  33. ** X Window graphics initially programmed 10/23-29/1991.
  34. ** PostScript graphics initially programmed 11/29-30/1992.
  35. ** Last code change made 9/22/1996.
  36. */
  37.  
  38. #include "astrolog.h"
  39.  
  40.  
  41. #ifdef MATRIX
  42. /*
  43. ******************************************************************************
  44. ** Assorted Calculations.
  45. ******************************************************************************
  46. */
  47.  
  48. /* Given a month, day, and year, convert it into a single Julian day value, */
  49. /* i.e. the number of days passed since a fixed reference date.             */
  50.  
  51. long MdyToJulian(mon, day, yea)
  52. int mon, day, yea;
  53. {
  54. #ifndef PLACALC
  55.   long im, j;
  56.  
  57.   im = 12*((long)yea+4800)+(long)mon-3;
  58.   j = (2*(im%12) + 7 + 365*im)/12;
  59.   j += (long)day + im/48 - 32083;
  60.   if (j > 2299171)                   /* Take care of dates in */
  61.     j += im/4800 - im/1200 + 38;     /* Gregorian calendar.   */
  62.   return j;
  63. #else
  64.   int fGreg = fTrue;
  65.  
  66.   if (yea < yeaJ2G || (yea == yeaJ2G &&
  67.     (mon < monJ2G || (mon == monJ2G && day < 15))))
  68.     fGreg = fFalse;
  69.   return (long)RFloor(julday(mon, day, yea, 12.0, fGreg)+rRound);
  70. #endif
  71. }
  72.  
  73.  
  74. /* Like above but return a fractional Julian time given the extra info. */
  75.  
  76. real MdytszToJulian(mon, day, yea, tim, dst, zon)
  77. int mon, day, yea;
  78. real tim, dst, zon;
  79. {
  80.   return (real)MdyToJulian(mon, day, yea) +
  81.     (DecToDeg(tim) + DecToDeg(zon) - DecToDeg(dst)) / 24.0;
  82. }
  83.  
  84.  
  85. /* Take a Julian day value, and convert it back into the corresponding */
  86. /* month, day, and year.                                               */
  87.  
  88. void JulianToMdy(JD, mon, day, yea)
  89. real JD;
  90. int *mon, *day, *yea;
  91. {
  92. #ifndef PLACALC
  93.   long L, N, IT, JT, K, IK;
  94.  
  95.   L  = (long)RFloor(JD+rRound)+68569L;
  96.   N  = Dvd(4L*L, 146097L);
  97.   L  -= Dvd(146097L*N + 3L, 4L);
  98.   IT = Dvd(4000L*(L+1L), 1461001L);
  99.   L  -= Dvd(1461L*IT, 4L) - 31L;
  100.   JT = Dvd(80L*L, 2447L);
  101.   K  = L-Dvd(2447L*JT, 80L);
  102.   L  = Dvd(JT, 11L);
  103.   JT += 2L - 12L*L;
  104.   IK = 100L*(N-49L) + IT + L;
  105.   *mon = (int)JT; *day = (int)K; *yea = (int)IK;
  106. #else
  107.   double tim;
  108.  
  109.   revjul(JD, JD >= 2299171.0 /* October 15, 1582 */, mon, day, yea, &tim);
  110. #endif
  111. }
  112.  
  113.  
  114. /* This is a subprocedure of CastChart(). Once we have the chart parameters, */
  115. /* calculate a few important things related to the date, i.e. the Greenwich  */
  116. /* time, the Julian day and fractional part of the day, the offset to the    */
  117. /* sidereal, and a couple of other things.                                   */
  118.  
  119. real ProcessInput(fDate)
  120. bool fDate;
  121. {
  122.   real Off, Ln;
  123.  
  124.   TT = RSgn(TT)*RFloor(RAbs(TT))+RFract(RAbs(TT))*100.0/60.0 +
  125.     (DecToDeg(ZZ) - DecToDeg(SS));
  126.   OO = DecToDeg(OO);
  127.   AA = Min(AA, 89.9999);        /* Make sure the chart isn't being cast */
  128.   AA = Max(AA, -89.9999);       /* on the precise north or south pole.  */
  129.   AA = RFromD(DecToDeg(AA));
  130.  
  131.   /* if parameter 'fDate' isn't set, then we can assume that the true time */
  132.   /* has already been determined (as in a -rm switch time midpoint chart). */
  133.  
  134.   if (fDate) {
  135.     is.JD = (real)MdyToJulian(MM, DD, YY);
  136.     if (!us.fProgress || us.fSolarArc)
  137.       is.T = (is.JD + TT/24.0 - 2415020.5) / 36525.0;
  138.     else {
  139.       /* Determine actual time that a progressed chart is to be cast for. */
  140.  
  141.       is.T = ((is.JD + TT/24.0 + (is.JDp - (is.JD + TT/24.0)) / us.rProgDay) -
  142.         2415020.5) / 36525.0;
  143.     }
  144.   }
  145.  
  146.   /* Compute angle that the ecliptic is inclined to the Celestial Equator */
  147.   is.OB = RFromD(23.452294-0.0130125*is.T);
  148.  
  149.   Ln = Mod((933060-6962911*is.T+7.5*is.T*is.T)/3600.0);  /* Mean lunar node */
  150.   Off = (259205536.0*is.T+2013816.0)/3600.0;             /* Mean Sun        */
  151.   Off = 17.23*RSin(RFromD(Ln))+1.27*RSin(RFromD(Off))-(5025.64+1.11*is.T)*is.T;
  152.   Off = (Off-84038.27)/3600.0;
  153.   is.rSid = (us.fSidereal ? Off : 0.0) + us.rZodiacOffset;
  154.   return Off;
  155. }
  156.  
  157.  
  158. /* Convert polar to rectangular coordinates. */
  159.  
  160. void PolToRec(A, R, X, Y)
  161. real A, R, *X, *Y;
  162. {
  163.   if (A == 0.0)
  164.     A = rSmall;
  165.   *X = R*RCos(A);
  166.   *Y = R*RSin(A);
  167. }
  168.  
  169.  
  170. /* Convert rectangular to polar coordinates. */
  171.  
  172. void RecToPol(X, Y, A, R)
  173. real X, Y, *A, *R;
  174. {
  175.   if (Y == 0.0)
  176.     Y = rSmall;
  177.   *R = RSqr(X*X + Y*Y);
  178.   *A = Angle(X, Y);
  179. }
  180.  
  181.  
  182. /* Convert rectangular to spherical coordinates. */
  183.  
  184. real RecToSph(B, L, O)
  185. real B, L, O;
  186. {
  187.   real R, Q, G, X, Y, A;
  188.  
  189.   A = B; R = 1.0;
  190.   PolToRec(A, R, &X, &Y);
  191.   Q = Y; R = X; A = L;
  192.   PolToRec(A, R, &X, &Y);
  193.   G = X; X = Y; Y = Q;
  194.   RecToPol(X, Y, &A, &R);
  195.   A += O;
  196.   PolToRec(A, R, &X, &Y);
  197.   Q = RAsin(Y);
  198.   Y = X; X = G;
  199.   RecToPol(X, Y, &A, &R);
  200.   if (A < 0.0)
  201.     A += 2*rPi;
  202.   G = A;
  203.   return G;  /* We only ever care about and return one of the coordinates. */
  204. }
  205.  
  206.  
  207. /* Do a coordinate transformation: Given a longitude and latitude value,    */
  208. /* return the new longitude and latitude values that the same location      */
  209. /* would have, were the equator tilted by a specified number of degrees.    */
  210. /* In other words, do a pole shift! This is used to convert among ecliptic, */
  211. /* equatorial, and local coordinates, each of which have zero declination   */
  212. /* in different planes. In other words, take into account the Earth's axis. */
  213.  
  214. void CoorXform(azi, alt, tilt)
  215. real *azi, *alt, tilt;
  216. {
  217.   real x, y, a1, l1;
  218.   real sinalt, cosalt, sinazi, sintilt, costilt;
  219.  
  220.   sinalt = RSin(*alt); cosalt = RCos(*alt); sinazi = RSin(*azi);
  221.   sintilt = RSin(tilt); costilt = RCos(tilt);
  222.  
  223.   x = cosalt * sinazi * costilt;
  224.   y = sinalt * sintilt;
  225.   x -= y;
  226.   a1 = cosalt;
  227.   y = cosalt * RCos(*azi);
  228.   l1 = Angle(y, x);
  229.   a1 = a1 * sinazi * sintilt + sinalt * costilt;
  230.   a1 = RAsin(a1);
  231.   *azi = l1; *alt = a1;
  232. }
  233.  
  234.  
  235. /* This is another subprocedure of CastChart(). Calculate a few variables */
  236. /* corresponding to the chart parameters that are used later on. The      */
  237. /* astrological vertex (object number nineteen) is also calculated here.  */
  238.  
  239. void ComputeVariables(vtx)
  240. real *vtx;
  241. {
  242.   real R, R2, B, L, O, G, X, Y, A;
  243.  
  244.   is.RA = RFromD(Mod((6.6460656+2400.0513*is.T+2.58E-5*is.T*is.T+TT)*15.0-OO));
  245.   R2 = is.RA; O = -is.OB; B = AA; A = R2; R = 1.0;
  246.   PolToRec(A, R, &X, &Y);
  247.   X *= RCos(O);
  248.   RecToPol(X, Y, &A, &R);
  249.   is.MC = Mod(is.rSid + DFromR(A));           /* Midheaven */
  250. #if FALSE
  251.   L = R2;
  252.   G = RecToSph(B, L, O);
  253.   is.Asc = Mod(is.rSid + Mod(G+rPiHalf));     /* Ascendant */
  254. #endif
  255.   L= R2+rPi; B = rPiHalf-RAbs(B);
  256.   if (AA < 0.0)
  257.     B = -B;
  258.   G = RecToSph(B, L, O);
  259.   *vtx = Mod(is.rSid + DFromR(G+rPiHalf));    /* Vertex */
  260. }
  261.  
  262.  
  263. /*
  264. ******************************************************************************
  265. ** House Cusp Calculations.
  266. ******************************************************************************
  267. */
  268.  
  269. /* The following three functions calculate the Midheaven, Ascendant, and  */
  270. /* East Point of the chart in question, based on time and location. The   */
  271. /* first two are also used in some of the house cusp calculation routines */
  272. /* as a quick way to get the 10th and 1st cusps. The East Point object is */
  273. /* technically defined as the Ascendant's position at zero latitude.      */
  274.  
  275. real CuspMidheaven()
  276. {
  277.   real MC;
  278.  
  279.   MC = RAtn(RTan(is.RA)/RCos(is.OB));
  280.   if (MC < 0.0)
  281.     MC += rPi;
  282.   if (is.RA > rPi)
  283.     MC += rPi;
  284.   return Mod(DFromR(MC)+is.rSid);
  285. }
  286.  
  287. real CuspAscendant()
  288. {
  289.   real Asc;
  290.  
  291.   Asc = Angle(-RSin(is.RA)*RCos(is.OB)-RTan(AA)*RSin(is.OB), RCos(is.RA));
  292.   return Mod(DFromR(Asc)+is.rSid);
  293. }
  294.  
  295. real CuspEastPoint()
  296. {
  297.   real EP;
  298.  
  299.   EP = Angle(-RSin(is.RA)*RCos(is.OB), RCos(is.RA));
  300.   return Mod(DFromR(EP)+is.rSid);
  301. }
  302.  
  303.  
  304. /* These are various different algorithms for calculating the house cusps: */
  305.  
  306. real CuspPlacidus(deg, FF, fNeg)
  307. real deg, FF;
  308. bool fNeg;
  309. {
  310.   real LO, R1, XS, X;
  311.   int i;
  312.  
  313.   R1 = is.RA+RFromD(deg);
  314.   X = fNeg ? 1.0 : -1.0;
  315.   /* Looping 10 times is arbitrary, but it's what other programs do. */
  316.   for (i = 1; i <= 10; i++) {
  317.  
  318.     /* This formula works except at 0 latitude (AA == 0.0). */
  319.  
  320.     XS = X*RSin(R1)*RTan(is.OB)*RTan(AA == 0.0 ? 0.0001 : AA);
  321.     XS = RAcos(XS);
  322.     if (XS < 0.0)
  323.       XS += rPi;
  324.     R1 = is.RA + (fNeg ? rPi-(XS/FF) : (XS/FF));
  325.   }
  326.   LO = RAtn(RTan(R1)/RCos(is.OB));
  327.   if (LO < 0.0)
  328.     LO += rPi;
  329.   if (RSin(R1) < 0.0)
  330.     LO += rPi;
  331.   return DFromR(LO);
  332. }
  333.  
  334. void HousePlacidus()
  335. {
  336.   int i;
  337.  
  338.   chouse[1] = Mod(is.Asc-is.rSid);
  339.   chouse[4] = Mod(is.MC+rDegHalf-is.rSid);
  340.   chouse[5] = CuspPlacidus(30.0, 3.0, fFalse) + rDegHalf;
  341.   chouse[6] = CuspPlacidus(60.0, 1.5, fFalse) + rDegHalf;
  342.   chouse[2] = CuspPlacidus(120.0, 1.5, fTrue);
  343.   chouse[3] = CuspPlacidus(150.0, 3.0, fTrue);
  344.   for (i = 1; i <= cSign; i++) {
  345.     if (i <= 6)
  346.       chouse[i] = Mod(chouse[i]+is.rSid);
  347.     else
  348.       chouse[i] = Mod(chouse[i-6]+rDegHalf);
  349.   }
  350. }
  351.  
  352. void HouseKoch()
  353. {
  354.   real A1, A2, A3, KN, D, X;
  355.   int i;
  356.  
  357.   A1 = RSin(is.RA)*RTan(AA)*RTan(is.OB);
  358.   A1 = RAsin(A1);
  359.   for (i = 1; i <= cSign; i++) {
  360.     D = Mod(60.0+30.0*(real)i);
  361.     A2 = D/rDegQuad-1.0; KN = 1.0;
  362.     if (D >= rDegHalf) {
  363.       KN = -1.0;
  364.       A2 = D/rDegQuad-3.0;
  365.     }
  366.     A3 = RFromD(Mod(DFromR(is.RA)+D+A2*DFromR(A1)));
  367.     X = Angle(RCos(A3)*RCos(is.OB)-KN*RTan(AA)*RSin(is.OB), RSin(A3));
  368.     chouse[i] = Mod(DFromR(X)+is.rSid);
  369.   }
  370. }
  371.  
  372. void HouseEqual()
  373. {
  374.   int i;
  375.  
  376.   for (i = 1; i <= cSign; i++)
  377.     chouse[i] = Mod(is.Asc-30.0+30.0*(real)i);
  378. }
  379.  
  380. void HouseCampanus()
  381. {
  382.   real KO, DN, X;
  383.   int i;
  384.  
  385.   for (i = 1; i <= cSign; i++) {
  386.     KO = RFromD(60.000001+30.0*(real)i);
  387.     DN = RAtn(RTan(KO)*RCos(AA));
  388.     if (DN < 0.0)
  389.       DN += rPi;
  390.     if (RSin(KO) < 0.0)
  391.       DN += rPi;
  392.     X = Angle(RCos(is.RA+DN)*RCos(is.OB)-RSin(DN)*RTan(AA)*RSin(is.OB),
  393.       RSin(is.RA+DN));
  394.     chouse[i] = Mod(DFromR(X)+is.rSid);
  395.   }
  396. }
  397.  
  398. void HouseMeridian()
  399. {
  400.   real D, X;
  401.   int i;
  402.  
  403.   for (i = 1; i <= cSign; i++) {
  404.     D = RFromD(60.0+30.0*(real)i);
  405.     X = Angle(RCos(is.RA+D)*RCos(is.OB), RSin(is.RA+D));
  406.     chouse[i] = Mod(DFromR(X)+is.rSid);
  407.   }
  408. }
  409.  
  410. void HouseRegiomontanus()
  411. {
  412.   real D, X;
  413.   int i;
  414.  
  415.   for (i = 1; i <= cSign; i++) {
  416.     D = RFromD(60.0+30.0*(real)i);
  417.     X = Angle(RCos(is.RA+D)*RCos(is.OB)-RSin(D)*RTan(AA)*RSin(is.OB),
  418.       RSin(is.RA+D));
  419.     chouse[i] = Mod(DFromR(X)+is.rSid);
  420.   }
  421. }
  422.  
  423. void HousePorphyry()
  424. {
  425.   real X, Y;
  426.   int i;
  427.  
  428.   X = is.Asc-is.MC;
  429.   if (X < 0.0)
  430.     X += rDegMax;
  431.   Y = X/3.0;
  432.   for (i = 1; i <= 2; i++)
  433.     chouse[i+4] = Mod(rDegHalf+is.MC+i*Y);
  434.   X = Mod(rDegHalf+is.MC)-is.Asc;
  435.   if (X < 0.0)
  436.     X += rDegMax;
  437.   chouse[1]=is.Asc;
  438.   Y = X/3.0;
  439.   for (i = 1; i <= 3; i++)
  440.     chouse[i+1] = Mod(is.Asc+i*Y);
  441.   for (i = 1; i <= 6; i++)
  442.     chouse[i+6] = Mod(chouse[i]+rDegHalf);
  443. }
  444.  
  445. void HouseMorinus()
  446. {
  447.   real D, X;
  448.   int i;
  449.  
  450.   for (i = 1; i <= cSign; i++) {
  451.     D = RFromD(60.0+30.0*(real)i);
  452.     X = Angle(RCos(is.RA+D), RSin(is.RA+D)*RCos(is.OB));
  453.     chouse[i] = Mod(DFromR(X)+is.rSid);
  454.   }
  455. }
  456.  
  457. real CuspTopocentric(deg)
  458. real deg;
  459. {
  460.   real OA, X, LO;
  461.  
  462.   OA = ModRad(is.RA+RFromD(deg));
  463.   X = RAtn(RTan(AA)/RCos(OA));
  464.   LO = RAtn(RCos(X)*RTan(OA)/RCos(X+is.OB));
  465.   if (LO < 0.0)
  466.     LO += rPi;
  467.   if (RSin(OA) < 0.0)
  468.     LO += rPi;
  469.   return LO;
  470. }
  471.  
  472. void HouseTopocentric()
  473. {
  474.   real TL, P1, P2, LT;
  475.   int i;
  476.  
  477.   chouse[4] = ModRad(RFromD(is.MC+rDegHalf-is.rSid));
  478.   TL = RTan(AA); P1 = RAtn(TL/3.0); P2 = RAtn(TL/1.5); LT = AA;
  479.   AA = P1; chouse[5] = CuspTopocentric(30.0) + rPi;
  480.   AA = P2; chouse[6] = CuspTopocentric(60.0) + rPi;
  481.   AA = LT; chouse[1] = CuspTopocentric(90.0);
  482.   AA = P2; chouse[2] = CuspTopocentric(120.0);
  483.   AA = P1; chouse[3] = CuspTopocentric(150.0);
  484.   AA = LT;
  485.   for (i = 1; i <= 6; i++) {
  486.     chouse[i] = Mod(DFromR(chouse[i])+is.rSid);
  487.     chouse[i+6] = Mod(chouse[i]+rDegHalf);
  488.   }
  489. }
  490.  
  491.  
  492. /*
  493. ******************************************************************************
  494. ** Planetary Position Calculations.
  495. ******************************************************************************
  496. */
  497.  
  498. /* Given three values, return them combined as the coefficients of a */
  499. /* quadratic equation as a function of the chart time.               */
  500.  
  501. real ReadThree(r0, r1, r2)
  502. real r0, r1, r2;
  503. {
  504.   return RFromD(r0 + r1*is.T + r2*is.T*is.T);
  505. }
  506.  
  507.  
  508. /* Another coordinate transformation. This is used by the ComputePlanets() */
  509. /* procedure to rotate rectangular coordinates by a certain amount.        */
  510.  
  511. void RecToSph2(AP, AN, IN, X, Y, G)
  512. real AP, AN, IN, *X, *Y, *G;
  513. {
  514.   real R, D, A;
  515.  
  516.   RecToPol(*X, *Y, &A, &R); A += AP; PolToRec(A, R, X, Y);
  517.   D = *X; *X = *Y; *Y = 0.0; RecToPol(*X, *Y, &A, &R);
  518.   A += IN; PolToRec(A, R, X, Y);
  519.   *G = *Y; *Y = *X; *X = D; RecToPol(*X, *Y, &A, &R); A += AN;
  520.   if (A < 0.0)
  521.     A += 2.0*rPi;
  522.   PolToRec(A, R, X, Y);
  523. }
  524.  
  525.  
  526. /* Calculate some harmonic delta error correction factors to add onto the */
  527. /* coordinates of Jupiter through Pluto, for better accuracy.             */
  528.  
  529. void ErrorCorrect(ind, x, y, z)
  530. int ind;
  531. real *x, *y, *z;
  532. {
  533.   real U, V, W, A, S0, T0[4], FPTR *pr;
  534.   int IK, IJ, irError;
  535.  
  536.   irError = rErrorCount[ind-oJup];
  537.   pr = (lpreal)&rErrorData[rErrorOffset[ind-oJup]];
  538.   for (IK = 1; IK <= 3; IK++) {
  539.     if (ind == oJup && IK == 3) {
  540.       T0[3] = 0.0;
  541.       break;
  542.     }
  543.     if (IK == 3)
  544.       irError--;
  545.     S0 = ReadThree(pr[0], pr[1], pr[2]); pr += 3;
  546.     A = 0.0;
  547.     for (IJ = 1; IJ <= irError; IJ++) {
  548.       U = *pr++; V = *pr++; W = *pr++;
  549.       A += RFromD(U)*RCos((V*is.T+W)*rPi/rDegHalf);
  550.     }
  551.     T0[IK] = DFromR(S0+A);
  552.   }
  553.   *x += T0[2]; *y += T0[1]; *z += T0[3];
  554. }
  555.  
  556.  
  557. /* Another subprocedure of the ComputePlanets() routine. Convert the final */
  558. /* rectangular coordinates of a planet to zodiac position and declination. */
  559.  
  560. void ProcessPlanet(ind, aber)
  561. int ind;
  562. real aber;
  563. {
  564.   real ang, rad;
  565.  
  566.   RecToPol(spacex[ind], spacey[ind], &ang, &rad);
  567.   planet[ind] = Mod(DFromR(ang) /*+ NU*/ - aber + is.rSid);
  568.   RecToPol(rad, spacez[ind], &ang, &rad);
  569.   if (us.objCenter == oSun && ind == oSun)
  570.     ang = 0.0;
  571.   ang = DFromR(ang);
  572.   while (ang > rDegQuad)    /* Ensure declination is from -90..+90 degrees. */
  573.     ang -= rDegHalf;
  574.   while (ang < -rDegQuad)
  575.     ang += rDegHalf;
  576.   planetalt[ind] = ang;
  577. }
  578.  
  579.  
  580. /* This is probably the heart of the whole program of Astrolog. Calculate  */
  581. /* the position of each body that orbits the Sun. A heliocentric chart is  */
  582. /* most natural; extra calculation is needed to have other central bodies. */
  583.  
  584. void ComputePlanets()
  585. {
  586.   real helioret[oNorm+1], heliox[oNorm+1], helioy[oNorm+1], helioz[oNorm+1];
  587.   real aber = 0.0, AU, E, EA, E1, M, XW, YW, AP, AN, IN, X, Y, G, XS, YS, ZS;
  588.   int ind = oSun, i;
  589.   OE *poe;
  590.  
  591.   while (ind <= (us.fUranian ? oNorm : cPlanet)) {
  592.     if (ignore[ind] && ind > oSun)
  593.       goto LNextPlanet;
  594.     poe = &rgoe[IoeFromObj(ind)];
  595.  
  596.     EA = M = ModRad(ReadThree(poe->ma0, poe->ma1, poe->ma2));
  597.     E = DFromR(ReadThree(poe->ec0, poe->ec1, poe->ec2));
  598.     for (i = 1; i <= 5; i++)
  599.       EA = M+E*RSin(EA);            /* Solve Kepler's equation */
  600.     AU = poe->sma;                  /* Semi-major axis         */
  601.     E1 = 0.01720209/(pow(AU, 1.5)*
  602.       (1.0-E*RCos(EA)));                     /* Begin velocity coordinates */
  603.     XW = -AU*E1*RSin(EA);                    /* Perifocal coordinates      */
  604.     YW = AU*E1*pow(1.0-E*E,0.5)*RCos(EA);
  605.     AP = ReadThree(poe->ap0, poe->ap1, poe->ap2);
  606.     AN = ReadThree(poe->an0, poe->an1, poe->an2);
  607.     IN = ReadThree(poe->in0, poe->in1, poe->in2); /* Calculate inclination  */
  608.     X = XW; Y = YW;
  609.     RecToSph2(AP, AN, IN, &X, &Y, &G);  /* Rotate velocity coords */
  610.     heliox[ind] = X; helioy[ind] = Y;
  611.     helioz[ind] = G;                    /* Helio ecliptic rectangtular */
  612.     X = AU*(RCos(EA)-E);                /* Perifocal coordinates for        */
  613.     Y = AU*RSin(EA)*pow(1.0-E*E,0.5);   /* rectangular position coordinates */
  614.     RecToSph2(AP, AN, IN, &X, &Y, &G);  /* Rotate for rectangular */
  615.     XS = X; YS = Y; ZS = G;             /* position coordinates   */
  616.     if (FBetween(ind, oJup, oPlu))
  617.       ErrorCorrect(ind, &XS, &YS, &ZS);
  618.     ret[ind] =                                        /* Helio daily motion */
  619.       (XS*helioy[ind]-YS*heliox[ind])/(XS*XS+YS*YS);
  620.     spacex[ind] = XS; spacey[ind] = YS; spacez[ind] = ZS;
  621.     ProcessPlanet(ind, 0.0);
  622. LNextPlanet:
  623.     ind += (ind == oSun ? 2 : (ind != cPlanet ? 1 : uranLo-cPlanet));
  624.   }
  625.  
  626.   spacex[oEar] = spacex[oSun];
  627.   spacey[oEar] = spacey[oSun];
  628.   spacez[oEar] = spacez[oSun];
  629.   planet[oEar] = planet[oSun]; planetalt[oEar] = planetalt[oSun];
  630.   ret[oEar] = ret[oSun];
  631.   heliox[oEar] = heliox[oSun]; helioy[oEar] = helioy[oSun];
  632.   helioret[oEar] = helioret[oSun];
  633.   spacex[oSun] = spacey[oSun] = spacez[oSun] =
  634.     planet[oSun] = planetalt[oSun] = heliox[oSun] = helioy[oSun] = 0.0;
  635.   if (us.objCenter == oSun) {
  636.     if (us.fVelocity)
  637.       for (i = 0; i <= oNorm; i++)    /* Use relative velocity */
  638.         ret[i] = RFromD(1.0);         /* if -v0 is in effect.  */
  639.     return;
  640.   }
  641. #endif /* MATRIX */
  642.  
  643.   /* A second loop is needed for geocentric charts or central bodies other */
  644.   /* than the Sun. For example, we can't find the position of Mercury in   */
  645.   /* relation to Pluto until we know the position of Pluto in relation to  */
  646.   /* the Sun, and since Mercury is calculated first, another pass needed.  */
  647.  
  648.   ind = us.objCenter;
  649.   for (i = 0; i <= oNorm; i++) {
  650.     helioret[i] = ret[i];
  651.     if (i != oMoo && i != ind) {
  652.       spacex[i] -= spacex[ind];
  653.       spacey[i] -= spacey[ind];
  654.       spacez[i] -= spacez[ind];
  655.     }
  656.   }
  657.   for (i = oEar; i <= (us.fUranian ? oNorm : cPlanet);
  658.     i += (i == oSun ? 2 : (i != cPlanet ? 1 : uranLo-cPlanet))) {
  659.     if ((ignore[i] && i > oSun) || i == ind)
  660.       continue;
  661.     XS = spacex[i]; YS = spacey[i]; ZS = spacez[i];
  662.     ret[i] = (XS*(helioy[i]-helioy[ind])-YS*(heliox[i]-heliox[ind]))/
  663.       (XS*XS+YS*YS);
  664.     if (ind == oEar)
  665.       aber = 0.0057756*RSqr(XS*XS+YS*YS+ZS*ZS)*DFromR(ret[i]); /* Aberration */
  666.     ProcessPlanet(i, aber);
  667.     if (us.fVelocity)                         /* Use relative velocity */
  668.       ret[i] = RFromD(ret[i]/helioret[i]);    /* if -v0 is in effect   */
  669.   }
  670.   spacex[ind] = spacey[ind] = spacez[ind] = 0.0;
  671. }
  672.  
  673.  
  674. #ifdef MATRIX
  675. /*
  676. ******************************************************************************
  677. ** Lunar Position Calculations
  678. ******************************************************************************
  679. */
  680.  
  681. /* Calculate the position and declination of the Moon, and the Moon's North  */
  682. /* Node. This has to be done separately from the other planets, because they */
  683. /* all orbit the Sun, while the Moon orbits the Earth.                       */
  684.  
  685. void ComputeLunar(moonlo, moonla, nodelo, nodela)
  686. real *moonlo, *moonla, *nodelo, *nodela;
  687. {
  688.   real LL, G, N, G1, D, L, ML, L1, MB, T1, Y, M = 3600.0, T2;
  689.  
  690.   T2 = is.T*is.T;
  691.   LL = 973563.0+1732564379.0*is.T-4.0*T2; /* Compute mean lunar longitude    */
  692.   G = 1012395.0+6189.0*is.T;              /* Sun's mean longitude of perigee */
  693.   N = 933060.0-6962911.0*is.T+7.5*T2;     /* Compute mean lunar node         */
  694.   G1 = 1203586.0+14648523.0*is.T-37.0*T2; /* Mean longitude of lunar perigee */
  695.   D = 1262655.0+1602961611.0*is.T-5.0*T2; /* Mean elongation of Moo from Sun */
  696.   L = (LL-G1)/M; L1 = ((LL-D)-G)/M;       /* Some auxiliary angles           */
  697.   T1 = (LL-N)/M; D = D/M; Y = 2.0*D;
  698.  
  699.   /* Compute Moon's perturbations. */
  700.  
  701.   ML = 22639.6*RSinD(L) - 4586.4*RSinD(L-Y) + 2369.9*RSinD(Y) +
  702.     769.0*RSinD(2.0*L) - 669.0*RSinD(L1) - 411.6*RSinD(2.0*T1) -
  703.     212.0*RSinD(2.0*L-Y) - 206.0*RSinD(L+L1-Y);
  704.   ML += 192.0*RSinD(L+Y) - 165.0*RSinD(L1-Y) + 148.0*RSinD(L-L1) -
  705.     125.0*RSinD(D) - 110.0*RSinD(L+L1) - 55.0*RSinD(2.0*T1-Y) -
  706.     45.0*RSinD(L+2.0*T1) + 40.0*RSinD(L-2.0*T1);
  707.  
  708.   *moonlo = G = Mod((LL+ML)/M+is.rSid);    /* Lunar longitude */
  709.  
  710.   /* Compute lunar latitude. */
  711.  
  712.   MB = 18461.5*RSinD(T1) + 1010.0*RSinD(L+T1) - 999.0*RSinD(T1-L) -
  713.     624.0*RSinD(T1-Y) + 199.0*RSinD(T1+Y-L) - 167.0*RSinD(L+T1-Y);
  714.   MB += 117.0*RSinD(T1+Y) + 62.0*RSinD(2.0*L+T1) -
  715.     33.0*RSinD(T1-Y-L) - 32.0*RSinD(T1-2.0*L) - 30.0*RSinD(L1+T1-Y);
  716.   *moonla = MB =
  717.     RSgn(MB)*((RAbs(MB)/M)/rDegMax-RFloor((RAbs(MB)/M)/rDegMax))*rDegMax;
  718.  
  719.   /* Compute position of the North Lunar Node, either True or Mean. */
  720.  
  721.   if (us.fTrueNode)
  722.     N = N+5392.0*RSinD(2.0*T1-Y)-541.0*RSinD(L1)-442.0*RSinD(Y)+
  723.       423.0*RSinD(2.0*T1)-291.0*RSinD(2.0*L-2.0*T1);
  724.   *nodelo = Mod(N/M+is.rSid);
  725.   *nodela = 0.0;
  726. }
  727. #endif /* MATRIX */
  728.  
  729. /* matrix.c */
  730.